function [spikes Gampa]=izhikevich(stim_conductance_exc, stim_conductance_inh)

C=100;vr=-60;vt=-40;k=3;
a=0.01;b=5;c=-50;d=100;
vpeak=35;

T=1000;

dt_ms=1;


spikes = 0;
n=round(T/dt_ms);
v=vr*ones(1,n); u=0*v;
si= stim_conductance_inh .*[zeros(1,20) ones(1,n-20)];
se= stim_conductance_exc .*[zeros(1,20) ones(1,n-20)];

% model conductances.
tau_ampa  = exp(-dt_ms/5);
tau_nmda  = exp(-dt_ms/150);
tau_gabaa = exp(-dt_ms/6);
tau_gabab = exp(-dt_ms/150);

E_AMPA	= 0.0;
E_NMDA	= 0.0;
E_GABAa	= -70;
E_GABAb = -90;
Ggabaa = 0;Ggabab = 0; Gampa = 0;Gnmda = 0;

Ggabaa=zeros(1,n);
Ggabab=zeros(1,n);
Gampa=zeros(1,n);
Gnmda=zeros(1,n);
I=zeros(1,n);

gain_ampa=1;
gain_nmda=1;
gain_gabaa=1;
gain_gabab=1;
% end model conductances.


for i = 1:n-1
    
    % conductances from synapses.
    Ggabaa(i+1) = Ggabaa(i)*tau_gabaa + dt_ms*(1-tau_gabaa)*gain_gabaa*si(i);
    Ggabab(i+1) = Ggabab(i)*tau_gabab + dt_ms*(1-tau_gabab)*gain_gabab*si(i);
    
    Gampa(i+1) = Gampa(i)*tau_ampa + dt_ms*(1-tau_ampa)*gain_ampa*se(i);
    Gnmda(i+1) = Gnmda(i)*tau_nmda + dt_ms*(1-tau_nmda)*gain_nmda*se(i);

    xx = (-80-v(i))/60;
	xx = xx.*xx;
	NMDAgate = xx./(1+xx);
    
     
    %g= ( Gampa(i+1)       +NMDAgate.*Gnmda(i+1)        + Ggabaa(i+1)         + Ggabab(i+1));
	%E= ( Gampa(i+1)*E_AMPA+NMDAgate.*Gnmda(i+1)*E_NMDA + Ggabaa(i+1)*E_GABAa + Ggabab(i+1)*E_GABAb);
    %I(i) = v(i).*g-E;
    I(i) = (v(i)-E_AMPA)*Gampa(i+1) + (v(i)-E_NMDA)*NMDAgate*Gnmda(i+1)+ ...
        (v(i)-E_GABAa)*Ggabaa(i+1)+(v(i)-E_GABAb)*Ggabab(i+1);
    
    
    v(i+1)=v(i)+.5*dt_ms*(k*(v(i)-vr)*(v(i)-vt)-u(i)-I(i))/C;
    v(i+1)=v(i+1)+.5*dt_ms*(k*(v(i+1)-vr)*(v(i+1)-vt)-u(i)-I(i))/C;
    u(i+1)=u(i)+dt_ms*a*(b*(v(i)-vr)-u(i));
    if v(i+1)>=vpeak
        v(i)=vpeak;
        v(i+1)=c;
        u(i+1)=u(i+1)+d;
        spikes = spikes + 1;
    end
end
figure(1)
subplot(3,1,1)
plot(dt_ms*(1:n),[v' u']);
subplot(3,1,2)
plot(dt_ms*(1:n),[Gampa' Gnmda' Ggabaa' Ggabab']);
title('conductances');
subplot(3,1,3);
plot(dt_ms*(1:n),I);
title('Isyn (pA)');
drawnow
